Repository navigation
Filter PacBio tumour-only SAVANA SVs against the SV PoN; cap ReConPlot SVs (#225) - #228
Open
ljwharbers wants to merge 6 commits into
Open
ljwharbers wants to merge 6 commits into
ljwharbers wants to merge 6 commits into
Conversation
…lot SVs SAVANA's PacBio classifier (`savana to --pb`) only rejects a call on matched-normal evidence, so without a normal nearly every germline SV comes out somatic, and RECONPLOT_SAVANA (one ggplot layer per SV) took 6-8 h (#225). The same VCF feeds VEP_SAVANA and Padfoot. - SAVANA_PON_FILTER / bin/savana_pon_filter.py: Severus's own PoN matching rule (extract_pon/add_pon) against the genome's pon_file, BND mates kept together. TUMORONLY_SAVANA routes PacBio somatic VCFs through it; ONT (SAVANA's TO model) and genomes without a PoN pass unchanged. --savana_pb_pon_filter turns it off. On a simulated PacBio tumour-only C8 (hg38) it keeps 761 of 27,591 events. - ReConPlot wrapper --max-svs (--reconplot_max_svs, default 2500): draw the best-supported junctions only; the harmonised SV table keeps them all. - Tests: SAVANA_PON_FILTER module test, TUMORONLY_SAVANA routing stub test, a RECONPLOT --max-svs case; RECONPLOT versions snapshots follow the VERSION bump. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Real-module run on the simulated PacBio tumour-only C8 (hg38): the PoN filter keeps 761 of 27,591 events and RECONPLOT_SAVANA finishes in 21 min with every figure drawn (unfiltered: not finished after 2 h). At --max-svs 2500 the genome-wide panel overflowed R's node stack, as it already does for BL1's 2,084 SAVANA junctions, so the default drops to 1000. The 761 survivors are still germline by C8's matched normal; the docs now say so. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
ljwharbers
marked this pull request as ready for review
November 2, 2026 12:56
![claude[bot]](https://avatars.githubusercontent.com/in/1236702?s=60&v=4){kind=small}
|
|
||
| ## after the tables, which keep every junction | ||
| n_sv <- nrow(sv) | ||
| sv <- cap_svs(sv, args$max_svs) |
There was a problem hiding this comment.
The cap is applied once across the whole genome, before regions are chosen, not "per figure".
cap_svs() runs on the genome-wide sv table before step 3 builds region_sets. Every run then draws from that same top-N set: per-chromosome (--layout separate), genome-wide, and the --regions/--genes focus panel.
- Per-chromosome figures: on a 27k-SV sample, each chromosome only gets its share of the genome-wide top 1000. A chromosome with 300 SVs, well under any overflow, may show only a few. Its title still says
top 1000 of 27591 SVs, counts that don't describe that figure. - Focus panel: a user-requested region whose SVs have low support can lose all of them, even if it has only a handful.
This contradicts nextflow.config ("draw at most this many SVs per figure") and the schema description ("per ReConPlot figure"). REVIEW.md requires help text to match actual behaviour.
The overflow only hits the genome-wide panel, so the fix should be one of:
- apply the cap per region set inside
plot_region_sets, after subsetting SVs to the region, and compute the title note there; or - apply it only to the
together/genome-wide run, and reword the docs to match.
| "type": "integer", | ||
| "description": "Mapping quality threshold for reads counted at a het-SNP site (SAVANA's --allele_mapq, default 5). Intended for test profiles running on downsampled data below that coverage floor." | ||
| }, | ||
| "savana_pb_pon_filter": { |
There was a problem hiding this comment.
docs/usage.md is missing both new parameters.
--savana_pb_pon_filter and --reconplot_max_svs are added to the schema, but not to docs/usage.md. REVIEW.md requires new parameters in both places.
The ReConPlot parameter table at docs/usage.md:466 already lists --reconplot_min_svlen, so --reconplot_max_svs belongs next to it. The SAVANA options section needs --savana_pb_pon_filter.
Comment on lines
+55
to
+59
| pon = {} | ||
| for key, entries in rows.items(): | ||
| entries.sort() | ||
| pon[key] = list(zip(*entries)) | ||
| return pon |
There was a problem hiding this comment.
If the PoN can't be parsed, the filter silently does nothing.
read_pon skips any line that has fewer than 7 comma-separated fields or non-integer columns. It does this without counting or warning. So if a genome's pon_file differs from the expected format, every row is skipped and pon comes back empty. Examples: tab-separated, or different column order or CI columns, like the CHM13 merged_PON_1000G_ASAP_withmedianCI.tsv.gz. The task then succeeds and publishes a pon_filtered VCF identical to the input. Downstream, VEP, Padfoot and ReConPlot get the unfiltered ~28k germline calls again.
Only the hand-made GRCh38-style fixture is tested. REVIEW.md asks for CHM13 to be checked as well.
Please fail, or at least warn on stderr, when no PoN rows were parsed, or when most data lines were skipped. Also confirm that the real hg38 and CHM13 PoN files parse to a non-zero row count.
|
Automated review: 3 findings (posted inline).
For a human reviewer to check by running things (this review could not run the pipeline or nf-test):
🤖 Generated with Claude Code |
A panel of normals in another layout (tab-separated, other columns) used to parse to nothing and pass every germline SV through as if filtered. Exit with an error when no row parses, warn when more than 1% of data lines are skipped, and log the parsed row count. Both shipped PoNs parse fully (hg38 1,039,498 rows; CHM13 1,490,071 plus its header); the C8 output is unchanged (761 of 27,591 events kept). Addresses review comment on PR #228. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The cap ran once on the genome-wide SV table before the region sets were built, so each per-chromosome figure got only its share of the genome-wide top N, a sparse focus region could lose all of its SVs, and every title carried the genome-wide count. Count the junctions each figure shows (either breakpoint in its regions) and cap only that figure; its title says top N of M with M for that figure. Uncapped figures still get the full table, so their output is unchanged. Addresses review comment on PR #228. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Addresses review comment on PR #228. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #225.
Why RECONPLOT_SAVANA took 6–8 h on PacBio tumour-only
SAVANA's PacBio classifier (
savana to --pb→classify_pacbio/pacbio_pass_somatic) only rejects a call on matched-normal evidence. With no normal, every germline SV with enough tumour reads is written to the somatic VCF. ReConPlot then adds one ggplot layer per SV. The same VCF also feeds VEP_SAVANA and Padfoot.Changes
SAVANA_PON_FILTER(bin/savana_pon_filter.py) filters the somatic VCF against the genome's Severus germline SV panel of normals (pon_file). It uses Severus's own matching rule (extract_pon/add_pon) and keeps or drops BND mates together.TUMORONLY_SAVANAsends only PacBio samples through it. ONT (SAVANA's trained tumour-only model) and genomes without a PoN pass through unchanged.<id>.pon_filtered.classified.somatic.vcfplus a summary TSV. The unfiltered VCF is still published.--savana_pb_pon_filter falseturns it off.--reconplot_max_svs(default 1000): ReConPlot draws only the best-supported SVs, and the title readstop N of M SVs. The harmonised*.reconplot_sv.tsvtable keeps every SV.Validation
Simulated PacBio tumour-only C8 (hg38): SAVANA's PacBio filter applied to the paired C8 breakpoints with the normal checks removed. The modules ran under Apptainer on one CPU.
--max-svs 1000--max-svs 2500devbecause of the same overflow.Tests
SAVANA_PON_FILTER.TUMORONLY_SAVANArouting: PacBio is filtered, ONT and no-PoN are not.--max-svs.versions.ymlmd5, from the wrapperVERSIONbump.PR checklist
docs/output.mdis updated.CHANGELOG.mdis updated.Follow-ups, not in this PR:
classify_pacbioupstream.🤖 Generated with Claude Code